Journal  of  Power  Sources  216  (2012)  152-161 


ELSEVIER 


Contents  lists  available  at  SciVerse  ScienceDirect 

Journal  of  Power  Sources 

journal  homepage:  www.elsevier.com/locate/jpowsour 


pri 

Sbb..«jtS 


A  quasi-three-dimensional  non-isothermal  dynamic  model  of  a  high-temperature 
proton  exchange  membrane  fuel  cell 

Jaeman  Park,  Kyoungdoug  Min* 

School  of  Mechanical  and  Aerospace  Engineering,  Seoul  National  University,  \  Gwanak-ro,  Gwanak-gu,  Seoul  151-744,  Republic  of  Korea 


HIGHLIGHTS 


►  The  model  can  predict  dynamic  characteristics  of  High-Temperature  PEM  Fuel  Cell. 

►  Increased  cell  temperature  results  in  enhanced  fuel  cell  performance. 

►  In  transient  state,  the  second  time  delay  is  related  to  the  MEA  temperature. 

►  The  transient  response  of  the  voltage  is  strongly  dependent  on  the  cell  temperature. 
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A  quasi-three-dimensional  non-isothermal  dynamic  model  of  a  high-temperature  proton  exchange 
membrane  fuel  cell,  operating  with  a  polybenzimidazole  membrane,  has  been  developed.  In  the  model, 
dynamic  equations  of  energy  conservation,  species  conservation  and  electrochemical  reaction  are  solved 
within  a  simplified  quasi-three-dimensional  geometry.  With  this  simplified  geometry,  the  model  can 
predict  the  distribution  of  fuel  cell  characteristics  and  physical  phenomena  in  a  high-temperature  proton 
exchange  membrane  fuel  cell  during  the  steady  and  the  transient  states  within  minutes,  without  using 
a  computational  fluid  dynamics  program.  The  simulation  results  of  the  current-voltage  polarization 
curve  are  validated  with  the  experimental  data  reported  in  the  literature.  The  simulation  results  show 
a  good  agreement  with  the  experimental  data  and  show  that  the  variation  of  temperature  strongly 
affects  the  fuel  cell’s  performance.  In  particular,  the  simulation  results  show  that  the  transient  response 
of  the  cell  voltage  is  dependent  upon  that  of  the  cell  temperature. 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

The  proton  exchange  membrane  fuel  cell  (PEMFC)  has  been 
regarded  as  a  promising  clean  energy  source  due  to  its  high  power 
density,  rapid  start-up,  high  efficiency  and  low  emissions  [1]. 
Although  great  advances  have  been  made  in  PEMFC  development, 
issues  such  as  water  management,  thermal  management,  sluggish 
electrochemical  cathode  kinetics  and  intolerance  to  impurities, 
such  as  CO,  remain  as  obstacles  to  commercialization  for  both 
transportation  and  stationary  applications  [2-4]. 

Recent  progress  in  PEMFC  technologies  has  focused  on  the 
need  to  develop  PEMFCs  that  operate  above  100  °C  [4-9].  A 
high-temperature  PEMFC  (HT-PEMFC)  with  a  polybenzimidazole 
(PBI)  membrane,  which  operates  in  the  temperature  range  of 
120—200  °C,  has  many  attractive  features  compared  with  typical 
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PEMFCs  [10,11].  One  of  the  advantages  of  the  PBI-based  HT-PEMFC 
is  that  such  a  system  does  not  require  humidification  [12,13].  The 
Nafion®  membrane,  which  is  based  on  perfluorosulfonated  poly¬ 
mers  and  is  used  in  typical  PEMFCs  with  an  operating  temperature 
range  of  50-80  °C,  requires  a  highly  hydrated  state  to  effectively 
conduct  protons.  Thus,  to  maintain  a  high  water  content  in  the  cell, 
external  humidification  is  required.  In  contrast,  a  PBI  membrane  is 
not  dependent  on  water  content  for  proton  conductivity  and  can 
operate  at  temperatures  up  to  200  °C;  therefore,  external  humidi¬ 
fication  is  not  required  and  the  product’s  water  exists  primarily 
in  the  vapor  phase  in  the  cell.  In  addition,  due  to  the  extremely  low 
electro-osmotic  drag  coefficient  of  the  PBI  membrane,  water 
transport  across  the  membrane  is  ignored  [14].  Consequently, 
water  management  is  no  longer  a  problem  in  the  HT-PEMFC  [15,16]. 
PBI-based  HT-PEMFCs  have  other  advantages  due  to  their  high 
operating  temperature,  including  improvement  in  the  electro¬ 
chemical  kinetics,  tolerance  of  the  poisonous  effects  of  CO  [17]  and 
the  application  of  a  simplified  cooling  system  [18,19]. 
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Nomenclature 

Re 

Reynolds  number  [-] 

R 

species  reaction  rate  (H2,  H20,  N2,  02)  [kmol  s-1] 

A 

surface  area  [m2] 

Rdif 

diffusion  resistance  [m3  s-1] 

C 

solid  specific  heat  capacity  [kj  kg  1  K-1] 

Sh 

Sherwood  number  [-] 

C 

species  molar  concentration  (H2>  H20,  N2,  02) 

t 

time  [s],  thickness  [m] 

Cv 

gas  constant- volume  specific  heat  capacity  [kj  kg-1  K 

T 

temperature  [I<] 

-1] 

V 

voltage  [V] 

Dm 

diffusion  coefficient  [m2  s  1] 

V 

fluid  velocity  [m  s-1] 

Dh 

hydraulic  diameter  [m] 

V 

volume  [m3] 

Ea 

activation  energy  [kj  kmol-1] 

X 

species  mole  fraction  (H2,  H20,  N2,  02)  [-] 

F 

Faraday’s  constant  [96,485  C  mol-1] 

f 

friction  factor  [-] 

Greek  letters: 

S  m 

mass  transport  coefficient  [m  s-1] 

0 

species  diffusion  flux  between  GDL  and  bulk  gases  (H2, 

g 

gravitational  acceleration  [ms-2] 

H20,  N2,  02)  [kmol  s-1] 

h 

enthalpy  [kj  kmol-1],  convective  heat  transfer 

a 

transfer  coefficient  [-] 

coefficient  [kW  m-2  K-1] 

£ 

porosity  of  GDL  [-] 

hf 

head  loss  [— ] 

K 

ionic  conductivity  [S  m-1] 

A  H 

enthalpy  of  formation  [kj  mol-1] 

P 

density  [kg  m-3] 

io 

exchange  current  density  [A  m-2] 

o 

electrode  conductivity  [S  m-1] 

I 

electrical  current  [A] 

kf 

fluid  conduction  heat  transfer  coefficient  [kW  m-1  K 

Superscripts  and  subscripts: 

-1] 

act 

activation  polarization 

ks 

solid  conduction  heat  transfer  coefficient  [kW  m  1  K 

GDL 

gas  diffusion  layer 

-1] 

g 

bulk  gas  control  volume 

L 

length  [m] 

in 

into  control  volume 

n 

number  of  participating  electrons  in  the  reaction  [-] 

MEA 

membrane  electrode  assembly 

N 

total  number  of  moles  [kmol] 

ohm 

ohmic  polarization 

N 

molar  flow  rate  [kmol  s-1] 

out 

out  of  control  volume 

Nud 

Nusselt  number  [-] 

ref 

reference  condition 

P 

pressure  [kPa] 

s 

solid  control  volume 

a 

heat  transfer  [kW] 

0 

standard  condition 

R 

fuel  cell  external/load  resistance  [Q] 

Even  though  the  HT-PEMFC  has  many  attractive  features,  only 
a  few  studies  on  HT-PEMFCs  have  been  performed.  In  particular, 
many  modeling  efforts  over  the  last  20  years  have  focused  on 
typical  PEMFCs  with  a  Nation®  membrane  [20-24];  however, 
research  on  HT-PEMFCs  with  a  PBI  membrane  began  approximately 
five  years  ago  and  there  remains  a  lack  of  modeling  work.  Only 
several  steady-state  models  for  the  HT-PEMFC  have  been  reported 
[14,25-31].  Cheddie  and  Munroe  [25-27]  were  the  first  authors 
to  develop  a  model  of  the  HT-PEMFC.  They  developed  the  one-, 
two-  and  three-dimensional  models.  These  models  could  predict 
the  oxygen  depletion  occur  in  the  cathode  catalyst  layer  and 
a  temperature  rise  of  up  to  20  K  at  a  power  density  of  1000  W  m  2. 
Peng  and  Lee  [14]  developed  a  three-dimensional  non-isothermal 
numerical  model  with  a  computational  fluid  dynamics  (CFD)  code 
and  they  found  that  thermal  effects  strongly  affect  the  fuel  cell 
performance.  Scott  et  al.  [28]  developed  a  steady-state  isothermal 
one-dimensional  model  which  considered  the  effect  of  water 
transport  across  the  membrane  on  the  conductivity  of  the  PBI 
membrane.  Ubong  et  al.  [29]  conducted  both  experimental  work 
and  numerical  simulations.  A  single  cell  with  an  active  area  of 
45  cm2  and  with  triple  serpentine  flow  channels  was  used  in  the 
experiment.  The  crossover  of  reactant  gases  through  the  membrane 
was  solved  analytically  by  Shamardina  et  al.  [30]  with  a  pseudo 
two-dimensional  steady-state  isothermal  model.  Lobato  et  al.  [31  ] 
developed  a  three-dimensional  model  that  was  solved  for  three 
different  flow  channel  geometries.  The  results  showed  that 
serpentine  flow  channels  achieved  the  best  performance.  In  addi¬ 
tion  to  these  steady-state  models,  a  few  dynamic  models  for 
HT-PEMFCs  were  published.  Zenith  et  al.  [33]  developed  a  dynamic 
model  using  the  characteristic  of  cell  resistance  as  an  input.  Peng 


et  al.  [34]  improved  upon  the  previous  model  [14]  with  an  accu¬ 
mulation  (unsteady)  term  and  observed  the  influence  of  the  charge 
double-layer  effect.  Sousa  et  al.  [35]  developed  a  dynamic 
two-dimensional  non-isothermal  model  and  simulated  the  tran¬ 
sient  response.  They  observed  the  double-layer  effect  and  degra¬ 
dation  effect  of  phosphoric  acid  loss  from  the  catalyst  layers  and 
platinum  sintering  in  the  cathode  catalyst  layer.  Additionally,  they 
reported  on  the  transient  response  of  the  oxygen  mole  fraction  and 
the  temperature  of  the  membrane  electrode  assembly  (MEA). 

The  models  that  have  been  developed  have  limitations.  First, 
the  active  area  of  the  MEA  is  small  (less  than  10  cm2).  With  this 
lab-scale  MEA  size,  previous  models  could  not  capture  the  local 
characteristics  of  the  HT-PEMFC.  Consideration  of  the  local  distri¬ 
bution  is  important  because  the  variations  in  the  local  character¬ 
istics  of  a  fuel  cell  may  cause  detrimental  effects  on  the 
performance  and  lifespan  of  the  fuel  cell.  Second,  most  of  the 
models  only  address  the  steady-state  characteristics.  To  understand 
the  transient  behavior  -  such  as  system  start-up,  shutdown  and 
rapid  changes  in  power  demand  -  a  dynamic  model  is  essential.  In 
fact,  the  models  based  on  CFD  may  not  be  manageable  or  suitable 
for  predicting  the  dynamic  behavior  and  developing  a  control 
system  because  they  require  large  computing  resources  and  have 
a  long  computational  time. 

The  objective  of  this  work  is  to  develop  a  quasi-three- 
dimensional  non-isothermal  dynamic  model  of  an  HT-PEMFC  with 
a  phosphoric  acid-doped  PBI  membrane.  This  work  is  the  only 
attempt  so  far  to  develop  a  model  of  an  HT-PEMFC  in  a  quasi-three- 
dimensional  geometry.  The  advantage  of  the  model  is  that  it  is 
simpler  in  geometry  compared  to  the  three-dimensional  models 
using  CFD  codes;  therefore,  it  is  possible  to  capture  local 
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distributions  of  temperature,  current  density  and  species  mole 
fractions  in  a  short  time.  The  calculation  time  was  123  s  for  the  model 
to  simulate  500  s  in  real  time  using  an  Intel®  Core™  Duo  Processor  at 
2.80  GHz  (Intel  Corporation,  Santa  Clara,  CA,  USA)  with  3.50  GB  RAM. 
Because  the  model  has  an  active  area  of  50  cm2,  it  is  possible  to 
consider  meaningful  simulation  results  of  local  distributions.  In 
addition,  this  model  can  be  used  to  capture  dynamic  responses 
during  a  step  change  in  current  and  to  elucidate  which  element 
could  affect  the  transient  response  by  considering  the  simulation 
results  from  the  transient  state.  Furthermore,  because  the  model  is 
implemented  using  MATLAB®-Simulink®  (MathWorks,  Inc.,  Natick, 
MA,  USA),  it  can  be  used  to  develop  and  test  control  strategies.  This 
model  could  also  be  expanded  to  include  integrated  fuel  cell 
systems.  Consequently,  this  model  can  be  useful  for  the  develop¬ 
ment  of  control  strategies  for  enhancing  overall  HT-PEMFC  system 
performance. 

2.  Model  description 

In  this  work,  a  dynamic  non-isothermal  quasi-three-dimen¬ 
sional  model  of  an  HT-PEMFC  with  a  PBI  membrane  is  developed. 
The  unit  cell  is  discretized  in  the  stream-wise  direction  (x,  y )  and  in 
the  cross-sectional  direction  (z).  By  this  method,  the  three- 
dimensional  geometry  of  the  model  can  be  simplified  (2+1 D, 
quasi-3D).  As  shown  in  Fig.  1,  the  unit  cell  is  discretized  along  the 
cross-sectional  direction  (z)  into  nine  control  volumes,  which 
consist  of  the  following:  (1)  the  solid  plate,  (2)  the  bulk  gas,  (3)  the 
GDL,  (4)  the  MEA  and  (5)  the  coolant.  These  control  volumes  are 
needed  to  understand  the  energy  and  the  mass  transport  of  local 
sections  of  the  cell  in  the  cross-sectional  direction.  In  addition  to 
this  discretization,  the  unit  cell  is  also  discretized  in  the  stream- 
wise  direction.  One  set  of  the  nine  control  volumes  in  the  cross- 
sectional  direction  is  repeated  in  the  stream-wise  direction.  This 
allows  the  model  to  capture  the  serpentine  flow  pattern,  and  the 
distributions  of  current  density,  temperature  and  species 


Flow  Inlet 


Fig.  1.  Illustration  of  the  stream-wise  and  the  cross-sectional  discretizations  of  an  HT- 
PEMFC 


Table  1 

Geometrical  HT-PEMFC  model  parameter  values. 


Description  geometry 

Value 

Units 

Cell  width  (x) 

0.07 

m 

Cell  height  (y) 

0.07 

m 

Depth  of  anode  gas  channel  (z) 

0.0015 

m 

Depth  of  cathode  gas  channel  (z) 

0.0015 

m 

Thickness  of  GDL  (z) 

0.0003 

m 

Thickness  of  electrolyte  (z) 

7  x  10-5 

m 

Thickness  of  separator  plates  (z) 

0.002 

m 

concentrations.  As  shown  in  Fig.  1,  the  model  contains  25  sets  of 
control  volumes  repeated  along  the  stream-wise  direction,  which 
represent  4-step  serpentine  flow  channels.  These  25  sets  of  control 
volumes  (a  total  of  225  control  volumes)  are  enough  to  capture  the 
local  characteristics  of  a  PEMFC  [20,21].  The  axial  heat  transfer 
network  is  considered  within  these  repeated  control  volumes. 
Tables  1  and  2  show  the  geometrical  and  general  model  parameter 
values,  respectively. 

2.1.  Assumptions 

The  major  assumptions  that  are  adopted  in  the  model  are  the 
same  as  that  of  Mueller  et  al.’s  work  [20].  In  HT-PEMFCs,  the  water 
in  the  cell  exists  only  in  the  vapor  phase,  so  a  single-phase  flow  is 
also  assumed  [27].  The  water  drag  coefficient  from  anode  to 
cathode  is  assumed  to  be  zero  because  of  the  property  of  the  PBI 
membrane  [29].  Additionally,  the  PBI  membrane  is  assumed  to  be 
impermeable  to  gas  flow  [27]. 

2.2.  Conservation  equations 

The  cell  characteristics  -  such  as  temperature,  species  mole 
fractions,  and  molar  flow  rate  of  each  control  volume  —  are  solved 
from  the  conservation  equations  [20,21  ].  In  addition,  in  the  present 
work,  some  equations  are  modified  to  adopt  the  characteristics  of 
HT-PEMFCs.  Energy  conservation  and  heat  transfer  equations  are 
used  to  solve  the  temperature,  and  the  species  conservation  and 
mass  transport  equations  are  used  to  solve  the  species  mole  frac¬ 
tions  and  the  molar  flow  rates.  These  equations  are  described  and 
given  in  the  following  sections  (2.2.1.,  2.2.2.). 


Table  2 

General  HT-PEMFC  model  parameter  values  for  the  validation  and  the  base  case. 


Description 

Value 

Units 

Transfer  coefficient 

1 

- 

GDL  porosity 

0.4 

- 

Cell  temperature 

423 

393  for  validation 

K 

Anode  stoichiometric  mass  flow 

1.5 

- 

rate  (H2) 

Cathode  stoichiometric  mass 

2.0 

- 

flow  rate  (air) 

Coolant  inlet  temperature 

423 

393  for  validation 

K 

Coolant  outlet  temperature 

425  @  0.7  A  cm-2 

394  @  0.4  A  cm-2  for 
validation 

K 

Anode/cathode  outlet  pressure 

101.325 

kPa 

Cathode  reference  exchange 

1.78  x  10“6 

Am"2 

current  density 

Tref,i 

353 

K 

Ei°lR 

16,456 

K 

Electrolyte  reference  conductivity 

1.87 

S  m  1 

Tref ,k 

423 

K 

E* /R 

4549.3 

K 

GDL  conductivity 

570 

S  m1 
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22.1.  Energy  conservation 

To  determine  the  temperature  within  the  control  volumes, 
energy  conservation  equations  are  solved.  The  temperatures  of 
each  solid  plate  control  volume  are  determined  by  solving  the 
following  ordinary  differential  equation  (ODE): 

^,vcS  =  E(i  in  (1) 

where  Qin  represents  heat  transfer  from  each  of  the  control 
volumes.  In  a  similar  way,  the  bulk  gas  and  the  coolant  temperature 
are  determined  by  solving  the  following  energy  equation: 

NCV-j^r  =  YAnhin  ~  Y,  ^out^out  +  Y,  Q-in  (2) 

Ninhin  and  Nou t/iout  in  Eq.  (2)  are  used  to  determine  the  change 
in  the  enthalpy  flux.  The  first  term  represents  the  enthalpy  flux  into 
the  control  volume  and  the  second  term  represents  the  enthalpy 
flux  out  of  the  control  volume. 

In  the  model,  the  temperature  of  the  control  volume  is  taken  to 
be  the  outlet  temperature  of  the  control  volume  and  this  outlet 
temperature  is  taken  to  be  the  inlet  temperature  of  the  next 
control  volume  in  the  stream-wise  direction.  Because  the  GDL  and 
the  MEA  are  assumed  to  have  a  lumped  temperature,  the 
temperature  of  the  GDL  and  MEA  is  found  by  combining  the  gas 
and  solid  control  volume  conservation  equations.  Furthermore, 
heat  generation  from  the  electrochemical  reactions  is  modeled  as 
follows: 


(E^a+E^cv)  g)§ 


E^in  -E^out h  out 
V-A  *  tj  I  V  I 

+  EQin+  HnF  1000 


(3) 


where  AH  is  the  enthalpy  of  the  formation  of  water,  and  V  and  I  are 
the  fuel  cell  voltage  and  current,  respectively. 


2.2.2.  Species  conservation 

The  species  molar  flow  rates  and  species  mole  fractions  at  each 
control  volume  exit  are  determined  from  species  mass  conserva¬ 
tion.  The  exit  mole  number  of  each  bulk  gas  control  volume  is  found 
from  the  species  conservation  equation: 

=  jvinXin  -  JVoutXout  +  ^  i  (4) 

The  exit  molar  flow  rate  is  determined  from  the  total  species 
conservation  equation: 

Nout  =  JVin  +  £5-£®  (5) 

O  in  Eqs.  (4)  and  (5)  is  the  species  diffusion  flux  from  the 
adjacent  GDL  control  volume.  The  “perfectly  stirred”  assumption 
[20]  is  applied  to  all  control  volumes  in  the  model  to  solve  the  exit 
mole  fraction  and  the  molar  flow  rate.  The  control  volume’s  mole 
fractions  are  assumed  to  be  those  at  the  control  volume  exits. 

The  species  mole  number  of  each  GDL  control  volume  is  deter¬ 
mined  as  follows: 

f=-E^  (6) 

where  R  is  the  local  electrochemical  reaction  rate  for  each  species. 
Because  there  is  no  water  transport  across  the  membrane,  the 
water  flux  from  the  MEA  is  ignored. 


2.3.  Heat  transfer 

Because  heat  is  generated  in  the  fuel  cell  due  to  irreversibility  in 
the  electrochemical  reactions,  the  extant  heat  transfer  between 
adjacent  control  volumes  EQ.in)  in  Eqs.  (1)— (3)  needs  to  be 
defined.  The  heat  generated  from  the  MEA  control  volume  is 
transferred  to  the  bulk  gas  through  convection  heat  transfer.  In 
addition,  some  amount  of  heat  from  the  MEA  control  volume  is  also 
transferred  to  the  solid  plate  through  conduction.  The  heat  in  the 
solid  plate  is  then  transferred  to  the  coolant  channel  (in  the  case  of 
a  cathode  solid  plate)  or  ambient  air  (in  the  case  of  an  anode  solid 
plate)  through  convection.  Furthermore,  some  amount  of  heat  from 
the  solid  plate  is  transferred  to  the  bulk  gas  because  the  solid  plate  is 
in  thermal  contact  with  the  bulk  gas.  In  addition  to  these  heat 
transfers  in  the  perpendicular  direction,  heat  is  transferred  between 
repeated  control  volumes  along  the  in-plane  direction.  Specifically, 
conduction  heat  transfer  between  adjacent  solid  plates  and  natural 
convection  at  the  edge  of  each  plate  are  included  in  the  model. 

Convection  heat  transfer  between  a  solid  and  a  gas  is  deter¬ 
mined  by  Newton’s  law  of  cooling: 

Q.  =  A  h(T2  -  7i)  (7) 

Conduction  heat  transfer  is  determined  by  Fourier’s  law: 

Q  =  Z  hi  (8) 


2.4.  Reaction  rates  and  species  diffusion 

Reaction  rates,  R,  and  species  diffusion,  &,  in  the  mass 
conservation  equations  (Eqs.  (4)-(6))  are  resolved  in  the  model. 

From  Faraday’s  law,  the  reaction  rates  of  both  half-reactions  are 
a  function  of  the  current  as  follows: 

^h2o  =  -%  =  — 2Rq2  =  jjfp  (9) 

Species  diffusion  is  considered  between  the  bulk  gas  and  GDL 
control  volumes.  The  driving  force  of  species  transport  from  the 
bulk  gas  to  the  GDL  is  convection  driven  by  a  concentration 
gradient  and  diffusion  in  the  GDL.  The  mass  transport  coefficient, 
gm,  at  the  interface  between  the  gas  channel  and  the  GDL  is 
obtained  based  on  the  Reynolds  analogy  between  heat  and  mass 
transfer: 


_  Sh  •  D  m 

= 


(10) 


where  Sh  is  the  Sherwood  number,  Dm  is  the  diffusion  coefficient 
and  Dh  is  the  hydraulic  diameter  of  the  gas  flow  channel.  The 
diffusion  coefficients  for  the  species  are  determined  by  the 
temperature  and  the  pressure.  A  modified  diffusion  coefficient, 
obtained  from  the  Bruggeman  correlation,  is  defined  to  account  for 
the  effects  of  porosity  and  tortuosity  in  the  GDL: 


Dm 


(11) 


Dm  =  a-5 -Dm  (12) 

where  D0  is  the  species  diffusion  coefficient  at  standard  pressure 

^eff 

and  temperature,  Dm  is  the  effective  species  diffusion  coefficient 
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and  e  is  the  GDL  porosity.  The  species  diffusion  flux  between  the 
bulk  gas  and  the  GDL  is  the  following: 

Rdif  =  A-  1  1  tGDL  (13) 

Sm  Dfft 

5  *  Rdif-(C2  -  C,)  (14) 

where  R  is  the  total  diffusion  resistance  of  each  species,  and  tGDL  is 
the  thickness  of  the  GDL 

2.5.  Electrochemical  model 


cell  length  was  considered.  The  Reynolds  number  was  calculated  to 
resolve  the  friction  factor  of  each  flow.  When  the  flow  had  a  Rey¬ 
nolds  number  under  2100,  it  was  assumed  to  be  laminar  flow 
where  the  Laminar  Friction  Constant,  /Re,  is  56.91  for  a  square  or 
a  rectangular  flow  geometry.  The  head  loss  can  be  calculated  by  Eq. 
(25)  and  the  pressure  drop  can  be  described  by  Eq.  (26).  The 
following  expressions  define  the  Reynolds  number,  friction  factor, 
head  loss  and  pressure  drop,  respectively  [32]: 


Re 


PvDh 

P 


(23) 


Re  <  2100  / 


56.91 

Re 


(24) 


The  operating  voltage  of  the  cell  is  determined  as  follows: 


V  —  VNernst;  ?7act  ^ohmic  (15) 

where  the  Nernst  voltage  over  100  °C  is  determined  based  on  the 
partial  pressures  in  the  GDL  and  on  the  MEA  temperature  [14].  The 
Nernst  voltage  is  found  by  the  following  equation: 


hf 


f 


DH2g 


Ap  =  pghf 


3.  Results  and  discussion 


(25) 

(26) 


^Nernst 


E°  +4.308  10  5  In 


pn2-(po  2)  °'5 


where 


(16) 


£°  =  1.17 -2.576- 10“4(r- 373.15)  (17) 

The  exchange  current  density  depends  on  the  temperature 
according  to  the  Arrhenius  equation  [31]: 


to  =  i'oefexp 


±(l _ M 

R  \T  T'ref.io ) 


(18) 


where  Ea°  is  the  activation  energy  and  igef  is  the  reference  exchange 
current  density  at  reference  temperature  T0.  The  activation  polari¬ 
zation  is  modeled  using  the  Tafel  equation: 


RT  , 

Vxt  —  «F  n 


/Wof\ 

y  to  c02  J 


(19) 


The  ionic  conductivity  of  the  PBI  membrane  is  described  by  the 
Arrhenius  equation  [31]: 


k  =  k0  exp 


e|/T 

R  \T 


(20) 


where  Ea  is  the  activation  energy  and  kq  is  the  pre-exponential 
factor. 

The  cell  resistance  is  modeled  by  the  membrane  ionic  conduc¬ 
tivity,  k ,  and  the  electrode  electric  conductivity,  a ,  and  the  ohmic 
polarization  is  proportional  to  the  cell  resistance: 


^ohmic  — 


fMEA  .  2tGDL 


(21) 


^ohmic  —  ^"^ohmic 


(22) 


3.1.  Model  validation 

The  polarization  curve  simulated  in  this  model  was  compared 
with  published  experimental  data  [31].  Because  the  model  in  this 
paper  was  developed  for  a  4-step  serpentine  flow  geometry,  the 
experimental  data  for  the  same  flow  geometry  were  selected.  Fuel 
cell  geometric  properties  and  conditions  of  the  experimental  work 
[31  ]  are  the  same  for  those  of  the  simulation  work  in  this  paper  as 
shown  in  Tables  1  and  2.  The  operating  temperature  for  the  vali¬ 
dation  is  393  K,  and  stoichiometric  mass  flow  ratio  of  anode  (FL) 
and  cathode  (air)  is  1.5  and  2.0,  respectively.  Since  the  current 
simulation  model  in  this  paper  is  non-isothermal  and,  yet,  there 
was  no  information  about  how  large  the  temperature  variation 
throughout  the  cell  in  the  experimental  work,  thus  it  was  assumed 
that  there  was  a  temperature  difference  of  1  I<  between  the  coolant 
inlet  and  the  coolant  outlet  at  the  load  current  condition  of 
0.4  A  cm-2.  Experimental  and  simulated  polarization  curves  over 
the  range  of  current  densities  from  0  to  0.4  A  cm-2  are  shown  in 
Fig.  2.  The  simulation  results  showed  a  good  agreement  with  the 
experimental  data  whose  error  range  is  within  4%. 

3.2.  Steady  results 

3.2.1.  Distribution  of  local  characteristics 

Because  the  model  is  intended  to  predict  local  characteristics 
such  as  current  density,  MEA  temperature  and  species  mole  frac¬ 
tion,  each  distribution  at  an  operating  temperature  of  423  K  is 
shown  in  Figs.  3-5.  These  results  were  simulated  at  a  high  current 
density  of  0.7  A  cm-2.  The  distribution  of  the  current  density  and 
the  MEA  temperature  were  investigated  with  two  different  coolant 
outlet  temperatures.  With  a  changing  coolant  flow  rate,  in  case  1, 
the  coolant  temperature  difference  between  the  inlet  and  the  outlet 
was  set  to  2  K  (423-425  K);  in  case  2,  the  coolant  temperature 
difference  was  set  to  7  K  (423-430  I<).  Simulation  results  are 
explained  below  and  summarized  in  Table  3. 


2.6.  Pressure  drop 

The  local  pressures  of  the  anode  and  the  cathode  channels  were 
calculated  to  predict  accurate  local  characteristics,  such  as  the 
Nernst  voltage.  The  pressure  loss  of  the  reactant  gas  flow  due  to  the 


32.1.1.  Species  mole  fraction.  The  local  hydrogen  mole  fraction  at 
the  anode  was  1  in  all  areas  because  the  cell  was  operating  with 
pure,  dry  hydrogen  and  no  water  diffusion  through  the  membrane. 
On  the  other  hand,  at  the  cathode,  the  oxygen  and  water  mole 
fractions  varied  along  the  flow  direction.  Fig.  3  shows  the  distri¬ 
butions  of  the  oxygen  mole  fraction  of  the  cathode  channel.  The 
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Fig.  2.  Comparison  of  polarization  curves  between  the  experimental  and  the  simula¬ 
tion  data  at  a  cell  temperature  of  393  K. 


oxygen  mole  fraction  at  the  cathode  channel  inlet  before  the  entry 
was  0.21  and  decreased  along  the  cathode  gas  channel.  The  water 
mole  fraction  at  the  cathode  channel  inlet  was  nearly  zero  because 
of  its  no-humidification  condition,  and  it  increased  along  the 
cathode  gas  channel  to  0.3.  The  tendency  of  the  variation  in  the 
water  mole  fraction  is  opposite  that  of  the  oxygen  mole  fraction 
(Fig.  3).  This  variation  is  because  oxygen  is  consumed  by  the  elec¬ 
trochemical  reaction  and  because  two  water  molecules  are 
produced  for  one  oxygen  molecule  consumed.  Thus,  the  water  mole 
fraction  increased  while  the  oxygen  mole  fraction  decreased  along 
the  cathode  gas  channel. 

3.2.12.  MEA  temperature.  The  MEA  temperature  increased  along 
the  coolant  channel  as  indicated  in  Fig.  4.  The  MEA  temperature 
increased  from  424.1  K  (at  the  cathode  inlet)  to  425.8  I<  (at  the 
cathode  outlet)  in  case  1  and  from  424.4  K  to  431.1  I<  in  case  2.  As 
heat  was  generated  by  the  cell,  the  heat  was  transferred  to  the 
coolant  and  the  coolant  temperature  increased.  Due  to  the 
increased  coolant  temperature,  the  amount  of  heat  removal  from 
the  MEA  to  the  coolant  decreased  along  the  coolant  channel.  This 
phenomenon  resulted  in  an  increase  in  the  MEA  temperature  along 
the  coolant  channel.  The  amount  of  heat  generated  by  the  cell  and 
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Fig.  3.  The  distribution  of  the  oxygen  mole  fraction  in  a  cathode  channel  at  a  cell 
operating  temperature  of  423  K  and  a  cell  current  density  of  0.7  A  cm-2. 


transferred  to  the  coolant  was  31.5  W  and  26.26  W,  respectively,  for 
the  heat  transfer  between  the  MEA  and  the  coolant  in  case  1.  In  case 
2,  the  amount  of  heat  generated  by  the  cell  and  transferred  to  the 
coolant  was  30.7  W  and  24.84  W,  respectively. 

32.1.3.  Current  density.  The  current  density  distributions  for  cases 
1  and  2  are  shown  in  Fig.  5.  In  case  1,  the  current  density  decreased 
along  the  cathode  gas  channel.  The  local  current  density  at  the 
cathode  outlet  (0.65  A  cm-2)  was  11%  lower  than  at  the  cathode 
inlet  (0.73  A  cm-2)  and  the  standard  deviation  of  the  current 
density  of  the  entire  area  was  0.025  A  cm-2.  This  result  was  caused 
by  a  variation  in  the  activation  overvoltage.  According  to  Eq.  (19), 
there  is  an  oxygen  concentration  term  and  the  activation  over¬ 
voltage  increased  while  the  oxygen  concentration  decreased  along 
the  cathode  gas  channel.  There  is  also  a  temperature  term,  such  as 
the  MEA  temperature  and  the  exchange  current  density,  in  the 
activation  overvoltage  equation.  Nevertheless,  as  described  in  the 
previous  section  (3.2.I.2.),  the  MEA  temperature  varied  only  0.4%  in 
case  1  and  this  change  did  not  significantly  affect  the  activation 
overvoltage.  For  the  same  reason,  the  membrane  ionic  conductivity 
varied  slightly  (1.922  S  m  1  at  the  cathode  inlet  and  2.008  S  m  1  at 
the  cathode  outlet).  Thus,  the  ohmic  resistance  also  varied  slightly 
and  this  effect  did  not  significantly  affect  the  local  current  density. 

Unlike  case  1,  in  case  2,  the  standard  deviation  of  the  current 
density  of  the  entire  area  was  0.009  A  cm-2  and  the  local  current 
density  was  quite  uniform  compared  to  case  1.  In  contrast  to  the 
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Fig.  5.  The  distribution  of  the  current  density  (A  cm  2)  at  a  cell  operating  temperature 
of  423  K  and  a  cell  current  density  of  0.7  A  cm-2,  a)  Case  1 ;  b)  Case  2. 


results  of  case  1,  the  temperature  influence  was  also  significant  for 
the  local  current  density  due  to  the  temperature  variance 
throughout  the  cell.  Specifically,  the  membrane  ionic  conductivity 
increased  along  the  cathode  gas  channel  (1.938  S  m-1  at  the 
cathode  inlet  and  2.287  S  m-1  at  the  cathode  outlet).  This  result  is 
because  the  ionic  conductivity  of  the  PBI  membrane  is  dependent 
on  temperature,  as  demonstrated  in  Eq.  (20),  and  the  ionic 
conductivity  is  proportional  to  the  MEA  temperature.  Thus,  the 
local  ohmic  resistance  in  the  cell,  which  is  strongly  dependent  upon 
the  membrane  ionic  conductivity,  decreases  as  the  MEA  tempera¬ 
ture  increases.  Consequently,  a  higher  temperature  improves  the 


Table  3 

Simulation  results  of  the  steady-state. 


Description 

Value 

Case  1 

Case  2 

Units 

Load  current 

0.7 

A  cm-2 

Operating  temperature 

423 

K 

Coolant  inlet  temperature 

423 

K 

Coolant  outlet  temperature 

425 

430 

K 

MEA  temperature  (cathode  inlet) 

424.1 

424.4 

K 

MEA  temperature  (cathode  outlet) 

425.  8 

431.1 

K 

Membrane  ionic  conductivity  (cathode  inlet) 

1.922 

1.938 

S  ITT1 

Membrane  ionic  conductivity  (cathode  outlet) 

2.008 

2.287 

S  ITT1 

Standard  deviation  of  current  density 

0.025 

0.009 

A  cm-2 

local  current  density.  Additionally,  in  the  activation  overvoltage,  the 
exchange  current  density  term  increases  as  the  MEA  temperature 
increases  as  in  Eq.  (18).  This  effect  could  reduce  the  activation 
overvoltage  along  the  cathode  gas  channel,  which  was  increased  by 
the  concentration  effect,  as  the  results  of  case  1  demonstrate.  Thus, 
increasing  temperature  enhances  fuel  cell  performance. 

3.3.  Dynamic  results 

3.3.1.  Cell  voltage  changes 

The  transient  response  of  the  voltage  when  a  step  change  in  the 
current  density  was  applied  to  the  model  is  shown  in  Fig.  6.  The 
current  density  was  instantaneously  increased  from  0.4  A  cm-2  to 
0.8  A  cm-2  at  500  s  and,  then,  instantaneously  decreased  from 
0.8  A  cnrr2  to  0.4  A  crrT2  at  700  s.  (The  mass  flow  rates  of  the 
reactant  gases  were  changed  immediately  by  a  corresponding 
amount  when  a  step  change  in  the  current  density  was  applied.) 
Undershoot  and  overshoot  occurred  in  the  voltage  in  response  to 
the  step  changes  in  the  current  density.  When  the  current  density 
was  changed  from  0.4  A  cm-2  to  0.8  A  cm-2,  the  voltage  rapidly 
decreased  from  0.48  V  to  0.3  V.  Then,  the  voltage  increased  to  its 
steady-state  value  of  0.33  V.  When  the  current  density  was  changed 
from  0.8  A  cm-2  to  0.4  A  cm-2,  the  voltage  suddenly  increased  from 
0.33  V  to  0.5  V  and,  then,  settled  to  its  steady-state  value  of  0.48  V. 
These  results  were  calculated  under  the  condition  at  which  coolant 
outlet  temperature  is  432  K  at  the  second  steady-state  (after  current 
step-up  from  0.4  A  cnrr2  to  0.8  A  cm-2)  while  the  coolant  inlet 
temperature  is  423  K  and  the  coolant  flow  rate  is  kept  constant 
despite  the  change  in  the  load.  Table  4  shows  simulation  conditions 
and  results  of  cell  voltage  changes. 

To  elucidate  the  voltage  transient  behavior,  especially  the 
undershoot  phenomenon,  the  transient  responses  of  the  average 
oxygen  mole  fraction  and  the  average  MEA  temperature  during  the 
current  step-up  were  considered.  The  transient  response  of  the 
average  oxygen  mole  fraction  in  the  catalyst  layer  during  the  load 
change  is  shown  in  Fig.  7(a).  When  the  load  was  changed  from 
0.4  A  cm-2  to  0.8  A  cm-2,  the  oxygen  mole  fraction  decreased 
smoothly  from  0.1119  to  0.0741  within  0.5  s.  These  results  can  be 
explained  by  the  following  phenomenon:  an  excess  of  oxygen  in 
the  catalyst  layer  was  consumed  immediately  after  the  load 
changed  and  the  change  in  the  local  oxygen  concentration  was 
delayed  because  the  oxygen  diffusion  cannot  follow  the  change  in 


Fig.  6.  Transient  response  of  voltage  as  the  current  density  changes  from  0.4  A  cm  2  to 
0.8  A  cm-2  at  500  s  and  from  0.8  A  cm-2  to  0.4  A  cm-2  at  700  s. 
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Table  4 

Simulation  results  of  the  transient  state. 


Description 

Value 

Units 

Current  step 

@  500  s 

0.4— 0.8 

A  cm-2 

@  700  s 

0.8— 0.4 

Coolant  inlet 

423 

K 

temperature 

Coolant  outlet 

@  1st  steady-state 

426 

K 

temperature 

@  2nd  steady-state 

432 

@  3rd  steady-state 

426 

Steady-state  value 

@  1st  steady-state 

0.48 

V 

@  2nd  steady-state 

0.33 

@  3rd  steady-state 

0.48 

current  density  instantaneously  [26].  This  phenomenon  resulted  in 
variations  in  the  open-circuit  voltage  and  the  activation  over¬ 
voltage.  As  described  in  Eqs.  (16)  and  (19),  the  open-circuit  voltage 
and  the  activation  overvoltage  depend  on  the  oxygen  mole  fraction 
at  the  cathode  catalyst  layer. 

The  transient  response  of  the  average  MEA  temperature,  when 
the  load  current  was  changed  from  0.4  A  cm-2  to  0.8  A  cm-2,  is 
shown  in  Fig.  7(b).  The  MEA  temperature  gradually  increased  from 
425.09  K  to  429.17  I<  within  25  s.  In  fuel  cells,  the  electrochemical 


reaction  is  more  active  in  the  high  current  density  area  and  the 
more  active  the  electrochemical  reaction  is,  the  more  heat  is 
released  from  the  cell.  Thus,  the  temperature  increased  during  the 
step  load  change.  Furthermore,  the  temperature  reached  the 
steady-state  value  after  a  few  tens  of  seconds  and  this  caused  a  time 
delay  in  reaching  the  steady  state  in  some  of  the  variables  that  are 
a  function  of  temperature.  Specifically,  the  activation  overvoltage 
and  the  ohmic  overvoltage,  which  directly  affect  the  I-V  curves  of 
the  fuel  cell,  varied  with  the  step  load  change. 

The  changes  of  the  activation  overvoltage  and  the  ohmic  over¬ 
voltage,  when  a  step  change  in  the  current  density  from  0.4  A  cm-2 
to  0.8  A  cm-2  occurred,  are  shown  in  Fig.  8.  The  open-circuit  voltage 
was  varied  by  only  0.04%  before  and  after  the  step  change  in  the 
current  density,  thus  it  is  not  considered  in  Fig.  8.  As  shown  in  Fig.  8, 
when  the  step  change  in  the  current  density  occurred,  the  activa¬ 
tion  overvoltage  and  the  ohmic  overvoltage  showed  somewhat 
different  responses  that  the  overvoltage  reached  their  highest  value 
during  the  overshoot  period.  The  activation  overvoltage,  which  is 
affected  by  both  the  oxygen  mole  fraction  and  the  MEA  tempera¬ 
ture,  reached  the  highest  value  of  0.568  V  at  0.3  s  after  the  step 
change  in  the  current  density  while  the  ohmic  overvoltage,  which 
is  affected  by  the  MEA  temperature,  reached  the  highest  value  of 
0.292  V  just  after  the  step  change  in  the  current  density  without 
any  delay.  After  that,  two  overvoltages  reached  their  steady-state 
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Fig.  7.  Transient  responses  of  the  oxygen  mole  fraction  and  the  MEA  temperature  at 
cathode  catalyst  layer  as  the  current  density  changes  from  0.4  A  cm-2  to  0.8  A  cm-2  at 
500  s  a)  oxygen  mole  fraction;  b)  MEA  temperature. 


Fig.  8.  Transient  responses  of  the  activation  overvoltage  and  the  ohmic  overvoltage  as 
the  current  density  changes  from  0.4  A  cm-2  to  0.8  A  cm-2  at  500  s  a)  activation 
overvoltage;  b)  ohmic  overvoltage. 
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value  within  a  few  tens  of  seconds.  These  behaviors  confirm  that 
the  effect  of  the  oxygen  mole  fraction  and  the  MEA  temperature  on 
the  overvoltage  responses  when  a  step  change  in  the  current 
density  occurred.  As  described  above,  the  transient  response  of 
oxygen  mole  fraction  affects  the  transient  response  of  the  activa¬ 
tion  overvoltage  in  the  early  stage  (within  a  few  tenths  of  seconds) 
and  the  transient  response  of  MEA  temperature  affects  the  activa¬ 
tion  overvoltage  and  the  ohmic  overvoltage  in  the  late  stage  (within 
a  few  tens  of  seconds).  Consequently,  the  transient  response  of 
voltage  is  affected  by  the  oxygen  mole  fraction  in  the  early  stage 
and  the  MEA  temperature  in  the  late  stage. 

The  tendencies  of  the  transient  responses  of  the  oxygen  mole 
fraction  and  the  MEA  temperature  were  also  reported  by  Sousa 
et  al.  [35].  They  reported  that  the  time  to  reach  steady-state  was 
within  0.15  s  for  the  oxygen  mole  fraction  and  200  s  for  the  MEA 
temperature.  Although  the  time  to  reach  the  steady-state  value  of 
the  MEA  temperature  is  somewhat  different  from  this  work,  it  is 
important  that  the  MEA  temperature  did  not  reach  the  steady- 
state  value  as  fast  as  the  oxygen  mole  fraction.  Additionally,  this 
tendency  was  compared  with  that  of  typical  PEMFCs  as  reported 
by  Cho  et  al.  [36].  They  demonstrated  that  there  are  two  different 
time  delays  in  the  transient  response  of  the  cell  voltage  when 
a  step  load  change  occurs.  The  first  time  delay  is  the  time  over 
which  the  voltage  decreases  until  it  reaches  the  minimum  level 
due  to  the  time  needed  for  oxygen  to  penetrate  the  GDL.  The 
second  time  delay  is  due  to  membrane  water  content  recovery 
which  is  related  to  the  membrane  ionic  conductivity.  In 
HT-PEMFCs,  there  were  also  two  different  time  delays  but  the 
reason  for  the  second  time  delay  differed  from  typical  PEMFCs.  The 
second  time  delay  in  HT-PEMFCs  is  due  to  the  time  needed  to  reach 
the  next  steady-state  value  of  the  MEA  temperature  which  is 
related  to  the  membrane  ionic  conductivity.  Consequently,  to 
maintain  performance  during  the  transient  state,  it  is  important  for 
an  HT-PEMFC  to  have  its  cell  temperature  reach  the  next  steady- 
state  value  quickly.  This  is  in  contrast  to  a  typical  PEMFC  in 
which  it  is  important  to  recover  membrane  water  content  quickly 
to  maintain  performance. 

3.3.2.  Temperature  influence 

The  transient  response  results  of  the  voltage  with  different 
coolant  outlet  temperatures  are  shown  in  Fig.  9.  There  are  three 
cases  that  control  the  degree  to  which  the  temperature  rises  in 


Time  /  s 

Fig.  9.  Transient  responses  of  voltage  with  different  coolant  outlet  temperatures  as  the 
current  density  changes  from  0.4  A  cm-2  to  0.8  A  cm-2  at  500  s. 


Table  5 

Simulation  results  of  the  transient  state  with  different  coolant  outlet  temperatures. 


Description 

Value 

Case  1 

Case  2 

Case  3 

Units 

Current  step 

0.4— 0.8 

A  cm-2 

Coolant  inlet  temperature 

423 

K 

Coolant  outlet  temperature 

432 

429 

426 

K 

Steady-state  value 

0.33 

0.31 

0.3 

V 

Time  to  reach  steady-state 

26 

17 

8 

s 

response  to  changes  in  the  coolant  flow  rate.  (Once  the  coolant  flow 
rate  was  determined  in  each  case,  it  did  not  change  despite  the 
change  in  the  load.)  Case  1  shows  the  greatest  difference,  9  K,  and 
case  2  and  case  3  show  differences  of  6  K,  and  3  K,  respectively, 
where  the  coolant  inlet  temperature  is  423  K.  Table  5  shows  these 
simulation  conditions.  When  the  current  density  was  changed  from 
0.4  A  cm-2  to  0.8  A  cm-2,  case  1  showed  the  best  steady-state 
performance  of  0.33  V;  however,  the  time  to  reach  the  steady 
state,  26  s,  was  the  longest  among  the  cases.  Case  2  showed 
a  steady-state  performance  of  0.31  V  and  a  time  delay  of  17  s.  Case  3 
showed  the  worst  steady-state  performance,  0.3  V;  however,  the 
time  to  reach  the  steady  state,  8  s,  was  the  shortest  among  the 
cases.  These  results  indicate  that  the  transient  response  of  the 
voltage  in  the  HT-PEMFC  is  strongly  dependent  on  cell  temperature. 
To  verify  this,  the  transient  response  of  the  average  MEA  temper¬ 
ature  to  the  coolant  outlet  temperature  was  plotted  and  is  shown  in 
Fig.  10.  When  the  step  load  change  occurred,  the  temperature 
gradually  increased  to  the  next  state  within  a  few  tens  of  seconds 
and  the  time  duration  increased  with  the  coolant  outlet  tempera¬ 
ture.  Furthermore,  as  the  temperature  reached  the  steady-state 
value  after  a  few  tens  of  seconds,  this  caused  a  time  delay  in 
reaching  the  steady  state  in  some  of  the  variables  that  are  a  func¬ 
tion  of  temperature.  The  activation  overvoltage  and  the  ohmic 
overvoltage,  which  directly  affect  the  I-V  curves  of  the  fuel  cell, 
varied  with  the  step  load  change.  Fig.  11  shows  the  variations  of  the 
activation  overvoltage  and  the  ohmic  overvoltage.  When  the  cell 
current  was  increased,  both  the  activation  overvoltage  and  the 
ohmic  overvoltage  increased  rapidly  and  decreased  to  the  next 
steady  state  over  a  few  tens  of  seconds.  These  results  indicate  why 
the  transient  response  of  the  voltage  varied  with  different  coolant 
outlet  temperatures. 


Fig.  10.  Transient  responses  of  the  MEA  temperature  as  the  current  density  changes 
from  0.4  A  cm-2  to  0.8  A  cm-2  at  500  s. 
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Fig.  11.  Transient  responses  of  the  activation  overvoltage  and  the  ohmic  overvoltage 
(lines  with  circles)  with  different  coolant  outlet  temperatures  as  the  current  density 
changes  from  0.4  A  cm-2  to  0.8  A  cm-2  at  500  s. 


4.  Conclusions 

A  dynamic  non-isothermal  quasi-three-dimensional  model  of 
a  high-temperature  PEMFC  with  a  PBI  membrane  has  been  devel¬ 
oped  and  simulated  at  various  operating  conditions  in  both  steady- 
state  and  transient  conditions. 

The  model  is  developed  using  the  MATLAB®-Simulink® 
programming  environment.  The  model  is  based  on  the  dynamic 
equations  of  mass  and  energy  conservation,  mass  and  heat  transfer 
and  electrochemical  reaction.  The  model  consists  of  control 
volumes  that  are  introduced  in  both  the  cross-sectional  direction 
and  the  stream-wise  direction.  With  these  control  volumes,  the 
model  can  predict  the  HT-PEMFC  characteristics.  Specifically,  the 
model  can  capture  distributions  of  local  current  density,  tempera¬ 
ture,  and  species  concentrations  within  a  few  minutes  without 
using  a  CFD  program.  To  validate  the  model,  the  polarization  curve 
of  the  simulation  data  is  compared  to  that  of  experimental  data 
from  the  literature  and  the  curve  shows  a  good  agreement  with  the 
experimental  data. 

The  simulation  results  of  the  steady-state  condition  showed  that 
there  is  a  strong  relationship  between  the  local  current  density  and 
the  local  MEA  temperature.  The  activation  overvoltage  increases 
due  to  oxygen  depletion  along  the  cathode  gas  channel  and  this 
causes  a  reduction  in  current  density.  The  rise  in  temperature 
through  the  cell  can  be  controlled  by  changing  the  coolant  flow  rate 
and  a  higher  temperature  enhances  membrane  ionic  conductivity. 
In  that  case,  the  ohmic  resistance  decreases  and  the  current  density 
increases.  These  trade-off  relationships  maintain  uniform  fuel  cell 
performance. 

The  dynamic  results  indicate  that  when  the  current  is  instan¬ 
taneously  increased  and  decreased,  a  voltage  undershoot  and 
overshoot  phenomenon  is  found.  This  undershoot  and  overshoot 
behavior  is  the  result  of  delays:  a  delay  in  the  local  oxygen 
concentration  because  the  oxygen  diffusion  cannot  follow  the 
change  in  current  density  instantaneously  and  a  delay  in  the  MEA 
temperature  reaching  the  steady-state  value.  The  time  delay  due  to 
oxygen  depletion  is  on  the  order  of  0.1  s  but  for  the  case  of  the  MEA 


temperature,  the  delay  time  is  on  the  order  of  10  s.  These  results 
indicate  that  the  transient  response  of  voltage  on  the  HT-PEMFC  is 
strongly  dependent  on  cell  temperature. 

Because  the  model  can  predict  dynamic  variations  of  charac¬ 
teristics  related  to  fuel  cell  performance  for  various  operating 
conditions  with  fast  simulation  times,  it  may  be  useful  for  the 
development  of  optimal  control  strategies  for  enhancing  HT-PEMFC 
performance.  In  the  future,  the  model  can  be  expanded  to  include 
the  integrated  fuel  cell  system. 
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